####################################################################################

#What You See and What You Get
#This file conducts the analysis for W2 survey
#Code to create Tables, Figures, and data cited in the main text and appendix 
#using Routputs created by running the previous three files (see documentation)


####################################################################################

rm(list=ls())
library(AER)
library(tidyverse)
library(readxl)
library(lubridate)
library(sandwich)
library(xtable)
library(corrgram)
library(forcats)
library(estimatr)
library(fastDummies)
library(randomizr)
library(ANOVAreplication)
library(miceadds)
library(clubSandwich)
library(texreg)
options(scipen=6)

source("Code/functions.R")

########################### Main Paper

#Main ITT and CACE effects
load("Routputs/out-W2-estimates-pooled.RData")

the.labels.plot_tp <-  c("PT Evaluation\nIndex", "Incumbent at Treat.\nEvaluation Index",
                         "Current Incumbent\nEvaluation Index", 
                         "Mobilization\nIndex",
                         "Attribution to\nLula, Dilma,\nor PT", 
                         "Attribute to\nPaes, Pezao,\nTemer, or MDB")

itts <- cbind(out.all$out.itt.lm$out.lmc 
              ,out.all$out.itt.lin$out.linc)  
itts <- itts[1:19,c(1,2,3,6,5,7,8,9,12,11,10)] 
itts_tp <- itts[1:6,] 
tp <- itts_tp

######################### Figure 1(b)

#Figure 1-b
png(file =  "Figures/fig-1-b.png",
    width = 8, height = 11, units = 'in', res = 600)
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
offsets <- c(0.12,-0.12)
ys <- -1:-nrow(tp)
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.7,.7)
par(mar=c(4,12,2,1))
plot(xs,c(-.75,-nrow(tp)-.25)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Intent-to-Treat Effects"
)
axis(side=2,at=ys
     ,labels=the.labels.plot_tp 
     ,las=2,tick=F,cex.axis=1.1, font = 2, 
     cex.axis = 1.3)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
abline(h=c(-6.5),col=gray(0.7))
#abline(h=c(-19.5),col=gray(0.7))
the.bgs <- c(1,"white")
the.cols <- c(1,1)
the.pchs <- c(21,21)
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i], cex = 1.1)
  points(tp[,estcol[i]],ys+offsets[i],bg=the.bgs[i],pch=the.pchs[i],cex=1.2,col=the.cols[i])	
}
legend(x="topright",legend=c("ITT+Ctrls (via Fixed Effects)","ITT+Ctrls (via Lin)")
       ,pch=the.pchs,col=the.cols,pt.bg=the.bgs,bty="n"
       ,cex=1.1,inset=c(-0.015,-0.015),xpd=NA)
dev.off()

##################### Online Appendix

################### Appendix C

load("Routputs/out-W2-balancebylottery.RData")
load("Routputs/out-W2-balancepooled.RData")

#Pooled balance
pol.bal <- xtable(bind_rows(pol.bal.stats$bal[c(1:4, 8)], 
                            pol.bal.stats$balr[c(1:4,8)]))

wald <- xtable(as.matrix(pol.bal.stats$wald, pol.bal.stats$pvalue))
pol.bal.stats$p.permutation #permutation p-value for Table C2
pol.bal.stats$pvalue #p-value for Table C2 (within table)

print(wald,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-c2.tex",sep=""))


print(pol.bal,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-c4.tex",sep=""))


#By lottery balance
ed0306 <- xtable(bind_rows(bal.stats$edital03.2011$balr[c(1:4, 8)], 
                           bal.stats$edital06.2011$balr[c(1:4,8)]))

wald03 <- xtable(bal.stats$edital03.2011$wald)
bal.stats$edital03.2011$p.permutation #permutation p-value for Table C6

wald06 <- xtable(bal.stats$edital06.2011$wald)
bal.stats$edital06.2011$p.permutation #permutation p-value for Table C6

print(ed0306,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-c9.tex",sep=""))

print(wald03,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-c6a.tex",sep=""))

print(wald06,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-c6b.tex",sep=""))


##################### Appendix D

#Objects
load("Routputs/out-W2-attrition-joint.RData")
load("Routputs/out-W2-attrition.RData")
load("Routputs/out-W2-benchmark.RData")

attritionW2 <- texreg(pol.attrition.stats, include.ci = F)

print(attritionW2,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-d4.tex",sep=""))


#Joint hypotheses test
attrition_overall_covi <- pol.attrition.joint.stats[[1]] 
attrition_answered_covi <- pol.attrition.joint.stats[[2]]
attrition_overall_cov <- pol.attrition.joint.stats[[3]]
attrition_answered_cov <- pol.attrition.joint.stats[[4]]


joint_overall <- xtable(waldtest(attrition_overall_covi, attrition_overall_cov))
joint_within <- xtable(waldtest(attrition_answered_covi, attrition_answered_cov))

print(joint_overall,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-d2a.tex",sep=""))

print(joint_within,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-d2b.tex",sep=""))

benchmarkW2 <- texreg(benchmarkW2, include.ci = F)

print(benchmarkW2,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tab-d5.tex",sep=""))

################## Appendix E

load("Routputs/out-W2-estimates-pooled.RData")

allc <- cbind(out.all$out.itt.lm$out.lmc 
              ,out.all$out.itt.lin$out.linc)  
allc <- allc[1:10,c(1,2,3,6,5,7,8,9,12,11,10)] 

the.labels.plot_tp <-  c("PT Eval. Index", "Incumbent at Treat.\nEval. Index",
                         "Current Incumbent\nEval Index", 
                         "Mobilization Index",
                         "Attribution to Lula/Dilma/PT", 
                         "Attribution to Paes/Pezao/Temer/MDB")

itts <- cbind(out.all$out.itt.lm$out.lmc 
            ,out.all$out.itt.lin$out.linc) 
itts_tp <- itts[1:6, c(1,2,3,6,5,7,8,9,12,11,10)] 
tp <- itts_tp

## Table for appendix (same results as the figure, which follows)
tpx <- tp
rownames(tpx) <- the.labels.plot_tp
tpx <- xtable(tpx,digits=3)
caption(tpx) <- "Analytical results from Figure 1b"
label(tpx) <- "tab:mainresults1b"
print(tpx, caption.placement = "top", 
            type="latex",
            hline.after=1,
            file="Tables/tab-e2.tex")

# Fig 1eb

itts <- cbind(out.all$out.itt.lin$out.lin
              ,out.all$out.itt.lin$out.linc)  
itts <- itts[1:19,c(1,2,3,6,5,7,8,9,12,11,10)] 
the.labels.plot_tp <-  c("PT Evaluation\nIndex", "Incumbent at Treat.\nEvaluation Index",
                         "Attribution to\nLula/Dilma/PT", 
                         "Attribution to Paes/\nPezao/Temer/MDB",
                         "Current Incumbent\nEvaluation Index", 
                         "Mobilization\nIndex")
itts_tp <- itts[c(1,2,5,6,3,4),]#original was [1:6]
tp <- itts_tp
pdf(file = paste("Figures/fig-e1b.pdf",sep=""),
    width = 8, height = 11)
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
offsets <- c(0.12,-0.12)
ys <- -1:-nrow(tp)
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.7,.7)
par(mar=c(4,12,2,1))
plot(xs,c(-.75,-nrow(tp)-.25)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Intent-to-Treat Effects"
)
axis(side=2,at=ys
     ,labels=the.labels.plot_tp 
     ,las=2,tick=F,cex.axis=1.1, font = 2)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
abline(h=c(-6.5),col=gray(0.7))
the.bgs <- c(1,"white")
the.cols <- c(1,1)
the.pchs <- c(21,21)
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i], cex = 1.1)
  points(tp[,estcol[i]],ys+offsets[i],bg=the.bgs[i],pch=the.pchs[i],cex=1.2,col=the.cols[i])	
}
legend(x="topright",legend=c("ITT  (via Lin)","ITT+Ctrls (via Lin)")
       ,pch=the.pchs,col=the.cols,pt.bg=the.bgs,bty="n"
       ,cex=1.1,inset=c(-0.015,-0.015),xpd=NA)
dev.off()

#### Figure E2

load("Routputs/out-W2-estimates-ipw.RData")

ittsipw <- rbind(out.all$out.itt.ipwc)
tp <- ittsipw

the.labels <-  c("PT Evaluation\nIndex", "Incumbent at Treat.\nEvaluation Index",
                 "Current Incumbent\nEvaluation Index", 
                 "Mobilization Index",
                 "Attribution to\nLula/Dilma/PT", 
                 "Attribution to\nPaes/Pezao/Temer/MDB",
                 "Lula Evaluation", "Dilma Evaluation",
                 "PT ID","Vote Haddad 2018",
                 "Paes Evaluation", 
                 "Bolsonaro Evaluation", "Crivela Evaluation", "Witzel Evaluation",
                 "Talk to Candidate","Support Candidate",
                 "Turnout 2016","Turnout 2018","Expected Turnout 2020", 
                 "Vote Pedro Paulo 16", "Vote Paes 2020", "Pezão Evaluation")

pdf("Figures/fig-e2a.pdf",
           width = 7, height = 10)
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
ys <- -1:-nrow(tp)
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.4,.4)
par(mar=c(4,10,2,1))
plot(xs,c(0,-nrow(tp))
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="ITTs"
     ,main="IPW Analysis")
axis(side=2,at=ys
     ,labels=the.labels 
     ,las=2,tick=F,cex.axis=0.8)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys,lty=3,col=gray(0.7))
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys)
  points(tp[,estcol[i]],ys,bg=0+i,pch=20+i,cex=0.8,col=0+i)	
}
legend(x="top",legend=c("ITT+Ctrls (via IPW)")
       ,pch=21,col=1,pt.bg=1,ncol=2,bty="n",
       cex=0.8,inset=-0.01,xpd=NA)
dev.off()

ittsipw <- rbind(out.all$out.itt.ipw)
tp <- ittsipw

pdf("Figures/fig-e2b.pdf",
           width = 7, height = 10)
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
ys <- -1:-nrow(tp)
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.4,.4)
par(mar=c(4,10,2,1))
plot(xs,c(0,-nrow(tp))
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="ITTs"
     ,main="IPW Analysis")
axis(side=2,at=ys
     ,labels=the.labels 
     ,las=2,tick=F,cex.axis=0.8)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys,lty=3,col=gray(0.7))
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys)
  points(tp[,estcol[i]],ys,bg=0+i,pch=20+i,cex=0.8,col=0+i)	
}
legend(x="top",legend=c("ITT (via IPW)")
       ,pch=21,col=1,pt.bg=1,ncol=2,bty="n",
       cex=0.8,inset=-0.01,xpd=NA)
dev.off()

####### Figures F
load("Routputs/out-W2-estimates-pooled.RData")

the.labels.plot_tp <-  c("Lula Evaluation", "Dilma Evaluation",
                         "PT ID","Vote Haddad 2018",
                         "Paes Evaluation", 
                         "Bolsonaro Evaluation", "Crivella Evaluation", "Witzel Evaluation",
                         "Talk to\nCandidate","Support\nCandidate",
                         "Turnout 2016","Turnout 2018","Expected \nTurnout 2020")

itts <- cbind(out.all$out.itt.lm$out.lmc 
              ,out.all$out.itt.lin$out.linc)  
itts_tp <- itts[7:19,c(1,2,3,6,5,7,8,9,12,11,10)] 
tp <- itts_tp
pdf(file = paste("Figures/fig-f1b.pdf",sep=""),
    width = 8, height = 11)
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
offsets <- c(0.12,-0.12)
ys <- -1:-nrow(tp)
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.7,.7)
par(mar=c(5,12,4,1))
plot(xs,c(-.85,-nrow(tp))
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Intent-to-Treat Effects"
     ,main="Early Lotteries Survey")
axis(side=2,at=ys
     ,labels=the.labels.plot_tp 
     ,las=2,tick=F,cex.axis=1.1, font = 2)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys-0.5,lty=3,col=gray(0.7))
the.bgs <- c(1,"white")
the.cols <- c(1,1)
the.pchs <- c(21,21)
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i], cex = 1.1)
  points(tp[,estcol[i]],ys+offsets[i],bg=the.bgs[i],pch=the.pchs[i],cex=1.2,col=the.cols[i])	
}
legend(x="topright",legend=c("ITT (via Fixed Effects)","ITT (via Lin)")
       ,pch=the.pchs,col=the.cols,pt.bg=the.bgs,bty="n"
       ,cex=1.1,inset=c(-0.015,-0.0375),xpd=NA)
dev.off()

#F2
the.labels.plot_tp <-  c("Lula Evaluation", "Dilma Evaluation",
                         "PT ID","Vote Haddad 2018",
                         "Paes Evaluation", 
                         "Bolsonaro Evaluation", "Crivella Evaluation", "Witzel Evaluation",
                         "Talk to\nCandidate","Support\nCandidate",
                         "Turnout 2016","Turnout 2018","Expected \nTurnout 2020")

itts <- cbind(out.all$out.itt.lm$out.lm 
              ,out.all$out.itt.lin$out.lin)  
itts_tp <- itts[7:19,c(1,2,3,6,5,7,8,9,12,11,10)] 
tp <- itts_tp
pdf(file = paste("Figures/fig-f2b.pdf",sep=""),
    width = 8, height = 11)
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
offsets <- c(0.12,-0.12)
ys <- -1:-nrow(tp)
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.7,.7)
par(mar=c(5,12,4,1))
plot(xs,c(-.85,-nrow(tp))
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Intent-to-Treat Effects"
     ,main="Early Lotteries Survey")
axis(side=2,at=ys
     ,labels=the.labels.plot_tp 
     ,las=2,tick=F,cex.axis=1.1, font = 2)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys-0.5,lty=3,col=gray(0.7))
the.bgs <- c(1,"white")
the.cols <- c(1,1)
the.pchs <- c(21,21)
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i], cex = 1.1)
  points(tp[,estcol[i]],ys+offsets[i],bg=the.bgs[i],pch=the.pchs[i],cex=1.2,col=the.cols[i])	
}
legend(x="topright",legend=c("ITT (via Fixed Effects)","ITT (via Lin)")
       ,pch=the.pchs,col=the.cols,pt.bg=the.bgs,bty="n"
       ,cex=1.1,inset=c(-0.015,-0.0375),xpd=NA)
dev.off()

##### Appendix G

the.labels.plot_tp <-  c("PT Evaluation\nIndex", 
                         "Incumbent at Treat.\nEvaluation Index", 
                         "Current Incumbent\nEvaluation Index", 
                         "Mobilization Index", 
                         "Attribution to\nLula/Dilma/PT", 
                         "Attribution to\nPaes/Pezão/Temer/MDB")

caces <- out.all$out.cace.lm$out.lmcacec 
caces <- caces[1:6,c(1,2,3,6,5)] 
tp <- caces
pdf(file = paste("Figures/fig-g1b.pdf",sep=""),
    width = 8, height = 11)
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
offsets <- c(0.12,-0.12)
ys <- -1:-nrow(tp)
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.8,.8)
par(mar=c(5,12,4,1))
plot(xs,c(-.75,-nrow(tp)-.25)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Complier Average Causal Effects"
)
axis(side=2,at=ys
     ,labels=the.labels.plot_tp 
     ,las=2,tick=F,cex.axis=1.1, font = 2)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
abline(h=c(-6.5),col=gray(0.7))
the.bgs <- c(1,"white")
the.cols <- c(1,1)
the.pchs <- c(21,21)
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i], cex = 1.1)
  points(tp[,estcol[i]],ys+offsets[i],bg=the.bgs[i],pch=the.pchs[i],cex=1.2,col=the.cols[i])	
}
dev.off()



